Workflow

Prepare Environment
Prepare Data
Disease Map
Point Pattern Analysis

4.1 Density based

4.2 Distance based

1. Prepare Environment

Required packages

# install.packages("sf")
# install.packages("tidyverse") #data wrangling
# install.packages("here") #working directory
# install.packages("janitor")
# install.packages("gtsummary")
# install.packages("DT")
# install.packages("stringr")
# install.packages("readxl")
# install.packages("broom")
# install.packages("mapview")
# install.packages("lubridate")
# install.packages("vctrs")
# install.packages("spatialECO") #for NNI
# install.packages("webshot2") #screenshot webpage

Load packages

library(sf)

## Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ purrr::%||%()   masks base::%||%()
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(here)

## here() starts at C:/Users/MY PC/OneDrive - Universiti Sains Malaysia/JKNK Spatial R course/spatial_workshop/drhazlienor.github.io

library(janitor)

## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test

library(gtsummary)
library(DT)
library(stringr)
library(readxl)
library(broom)
library(tmap)

## Breaking News: tmap 3.x is retiring. Please test v4, e.g. with
## remotes::install_github('r-tmap/tmap')

library(mapview)
library(lubridate)
library(maptools)

## Loading required package: sp

## Warning: multiple methods tables found for 'elide'

## Checking rgeos availability: TRUE
## Please note that 'maptools' will be retired during 2023,
## plan transition at your earliest convenience;
## some functionality will be moved to 'sp'.
## 
## Attaching package: 'maptools'
## 
## The following objects are masked from 'package:sp':
## 
##     elide, sp2Mondrian

library(spatstat)

## Loading required package: spatstat.data
## Loading required package: spatstat.univar
## spatstat.univar 3.0-0
## Loading required package: spatstat.geom
## spatstat.geom 3.3-2
## Loading required package: spatstat.random
## spatstat.random 3.3-1
## Loading required package: spatstat.explore
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## 
## The following object is masked from 'package:dplyr':
## 
##     collapse
## 
## spatstat.explore 3.3-1
## Loading required package: spatstat.model
## Loading required package: rpart
## spatstat.model 3.3-1
## Loading required package: spatstat.linnet
## spatstat.linnet 3.2-1
## 
## spatstat 3.1-1 
## For an introduction to spatstat, type 'beginner'

library(spdep)

## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

library(gridExtra)

## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

library(sparr)

## 
## 
## Welcome to
##    _____ ___  ____  ____  ____         
##   / ___// _ \/ _  \/ __ \/ __ \        
##   \__ \/ ___/ __  /  ___/  ___/        
##  ___/ / /  / / / / /\ \/ /\ \          
## /____/_/  /_/ /_/_/  \__/  \_\   v2.3-15
## 
## - type news(package="sparr") for an overview
## - type help("sparr") for documentation
## - type citation("sparr") for how to cite

library(grid)

## 
## Attaching package: 'grid'
## 
## The following object is masked from 'package:spatstat.geom':
## 
##     as.mask

library(spatialEco)

## 
## Attaching package: 'spatialEco'
## 
## The following object is masked from 'package:gridExtra':
## 
##     combine
## 
## The following objects are masked from 'package:spatstat.geom':
## 
##     is.empty, quadrats, shift
## 
## The following object is masked from 'package:spatstat.data':
## 
##     ants
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

library(plotly)

## 
## Attaching package: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout

library(gapminder)
library(webshot2)

2. Prepare Data

Types of data: polygon data, population data, point data
Import Data: Load spatial and attribute data (e.g., shapefiles, CSVs) into R using functions like read.csv(), sf::st_read() for vector data, or raster::raster() for raster data.
Clean Data: Handle missing values, correct data formats, and ensure the spatial data is properly projected using sf::st_transform() to standardize coordinate reference systems (CRS).
Explore Data: Use basic summaries (summary(), head()) and visual checks (plot(), View()) to understand data structure.

Polygon data

read polygon data - kelantan map

kel_map <- st_read(here("Map",
                    "kelantan.shp"))

## Reading layer `kelantan' from data source 
##   `C:\Users\MY PC\OneDrive - Universiti Sains Malaysia\JKNK Spatial R course\spatial_workshop\drhazlienor.github.io\Map\kelantan.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 66 features and 6 fields
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 371629.6 ymin: 503028.2 xmax: 519479.6 ymax: 690232.8
## Projected CRS: Kertau (RSO) / RSO Malaya (m)

st_geometry() extracts the geometric part (coordinates) of an sf object, allowing you to view or manipulate the spatial features

st_geometry(kel_map)

## Geometry set for 66 features 
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: 371629.6 ymin: 503028.2 xmax: 519479.6 ymax: 690232.8
## Projected CRS: Kertau (RSO) / RSO Malaya (m)
## First 5 geometries:

## POLYGON ((485501.8 669698.8, 485717.3 669694.6,...

## POLYGON ((487716.5 665649.5, 487615.4 665445.1,...

## POLYGON ((482744.8 660223.4, 482823.6 660137.8,...

## POLYGON ((486936.9 677358.5, 486990.5 677333.7,...

## POLYGON ((490841.7 668783.4, 490906.1 668691, 4...

Population data

read population data per subdistrict/mukim per year. You can get the population data from DOSM.

kel_mukim <- read_xlsx(here ("pop_kel.xlsx"))
kel_mukim %>% datatable()

Point data

read linelisting cointaining point data in .xlsx format

kel_lepto <- read_xlsx("leptospirosis.xlsx") %>% clean_names()
glimpse(kel_lepto)

## Rows: 1,106
## Columns: 26
## $ diagnosis                                  <chr> "Leptospirosis", "Leptospir…
## $ notifikasi_no                              <dbl> 2685725, 2728504, 2739963, …
## $ tahun_daftar                               <dbl> 2016, 2016, 2016, 2016, 201…
## $ epid_daftar                                <dbl> 6, 9, 11, 12, 18, 18, 19, 1…
## $ age                                        <dbl> 30, 23, 39, 43, 31, 34, 48,…
## $ alamat_semasa_kejadian                     <chr> "FELDA ARING", "LADANG U&I …
## $ poskod                                     <dbl> 18300, 18300, 18300, 18300,…
## $ latitude_rso                               <dbl> 478031, 459494, 441802, 488…
## $ longitude_rso                              <dbl> 548141, 564966, 547551, 547…
## $ latitude_wgs                               <dbl> 4.944824, 5.034372, 4.89743…
## $ longitude_wgs                              <dbl> 102.3491, 102.1477, 101.960…
## $ notifikasi_status                          <chr> "Daftar Kes", "Daftar Kes",…
## $ race                                       <chr> "Foreigner", "Foreigner", "…
## $ kewarganegaraan                            <chr> "Bukan Warganegara", "Bukan…
## $ gender                                     <chr> "Male", "Male", "Male", "Ma…
## $ nationality                                <chr> "INDONESIA", "INDONESIA", "…
## $ klasifikasi_kes                            <chr> "Sporadic", "Sporadic", "Sp…
## $ cara_pengesanan_kes                        <chr> "Pasif", "Pasif", "Pasif", …
## $ jenis_rawatan                              <chr> "Wad Perubatan", "Jabatan K…
## $ daerah                                     <chr> "GUA MUSANG", "GUA MUSANG",…
## $ mukim                                      <chr> "CHIKU", "CHIKU", "GALAS", …
## $ lokaliti                                   <chr> NA, "LADANG U & I, CIKU", "…
## $ diagnosis2                                 <chr> "LEPTOSPIROSIS", "LEPTOSPIR…
## $ sub_diagnosis                              <lgl> NA, NA, NA, NA, NA, NA, NA,…
## $ negeri                                     <chr> "KELANTAN", "KELANTAN", "KE…
## $ jenis_import_jangkitan_dalam_negara_negeri <chr> NA, NA, NA, NA, NA, NA, NA,…

convert all character to factor (categorical) to allow better understanding of the data

kel_lepto <- kel_lepto %>%
  mutate(across(where(is.character), as.factor))

explore the data

summary(kel_lepto)

##          diagnosis    notifikasi_no      tahun_daftar   epid_daftar   
##  Leptospirosis:1106   Min.   :2662339   Min.   :2016   Min.   : 1.00  
##                       1st Qu.:3005967   1st Qu.:2016   1st Qu.:14.00  
##                       Median :3510632   Median :2018   Median :28.00  
##                       Mean   :4701261   Mean   :2018   Mean   :27.57  
##                       3rd Qu.:5042478   3rd Qu.:2020   3rd Qu.:42.00  
##                       Max.   :9573673   Max.   :2022   Max.   :53.00  
##                                                                       
##       age                alamat_semasa_kejadian     poskod     
##  Min.   : 0.42   KG MUKA BUKIT      :   3       Min.   :11510  
##  1st Qu.:15.25   KG PANGLIMA BAYU   :   3       1st Qu.:16800  
##  Median :27.00   RPT SG PAS         :   3       Median :18000  
##  Mean   :30.95   FELCRA TERATAK BATU:   2       Mean   :17543  
##  3rd Qu.:45.00   FELDA ARING        :   2       3rd Qu.:18300  
##  Max.   :95.00   KAMPUNG BENDANG    :   2       Max.   :40200  
##                  (Other)            :1091       NA's   :151    
##   latitude_rso    longitude_rso     latitude_wgs   longitude_wgs  
##  Min.   :     6   Min.   : 63494   Min.   :4.611   Min.   :101.4  
##  1st Qu.:450171   1st Qu.:601283   1st Qu.:5.374   1st Qu.:102.1  
##  Median :464957   Median :638018   Median :5.764   Median :102.2  
##  Mean   :459729   Mean   :627313   Mean   :5.651   Mean   :102.1  
##  3rd Qu.:472980   3rd Qu.:665421   3rd Qu.:6.013   3rd Qu.:102.3  
##  Max.   :502469   Max.   :688574   Max.   :6.230   Max.   :102.5  
##  NA's   :736      NA's   :737                                     
##        notifikasi_status         race               kewarganegaraan
##  Abai Notifikasi:   5    Aborigines:  27   Bukan Warganegara:  45  
##  Daftar Kes     :1101    Chinese   :   9   Warganegara      :1061  
##                          Foreigner :  45                           
##                          Indian    :   5                           
##                          Malay     :1016                           
##                          Others    :   4                           
##                                                                    
##     gender        nationality   klasifikasi_kes cara_pengesanan_kes
##  Female:367   MALAYSIA  :1061   Outbreak:   5   Aktif:   6         
##  Male  :739   INDONESIA :  17   Sporadic:1101   Pasif:1100         
##               MYANMAR   :   9                                      
##               THAILAND  :   7                                      
##               BANGLADESH:   6                                      
##               INDIA     :   3                                      
##               (Other)   :   3                                      
##                         jenis_rawatan         daerah                mukim    
##  ICU                           : 58   GUA MUSANG :220   BATU MENGKEBANG: 89  
##  Jabatan Kecemasan & Kemalangan:385   KUALA KRAI :173   BERTAM         : 89  
##  Jabatan Pesakit Luar          :117   MACHANG    :143   GALAS          : 79  
##  Lain-lain                     : 23   PASIR MAS  :130   OLAK JERAM     : 67  
##  Wad Kanak-kanak               :122   KOTA BHARU :121   ULU SAT        : 62  
##  Wad Perubatan                 :401   TANAH MERAH: 99   CHIKU          : 52  
##                                       (Other)    :220   (Other)        :668  
##         lokaliti            diagnosis2   sub_diagnosis       negeri    
##  LATA REK   :   7   LEPTOSPIROSIS:1106   Mode:logical   KELANTAN:1106  
##  PERIA      :   7                        NA's:1106                     
##  TANAH PUTIH:   7                                                      
##  RPT SG PAS :   6                                                      
##  BATU BALAI :   5                                                      
##  (Other)    :1004                                                      
##  NA's       :  70                                                      
##  jenis_import_jangkitan_dalam_negara_negeri
##  KELANTAN:   3                             
##  NA's    :1103                             
##                                            
##                                            
##                                            
##                                            
##

identify observation with missing coordinates. you can update the missing coordinates in the same linelisting and re-run the command

# Identify rows with missing coordinate data
no_coordinates <- kel_lepto %>%
  filter(is.na(latitude_wgs) | is.na(longitude_wgs))

# Print rows with no coordinates
print(no_coordinates)

## # A tibble: 0 × 26
## # ℹ 26 variables: diagnosis <fct>, notifikasi_no <dbl>, tahun_daftar <dbl>,
## #   epid_daftar <dbl>, age <dbl>, alamat_semasa_kejadian <fct>, poskod <dbl>,
## #   latitude_rso <dbl>, longitude_rso <dbl>, latitude_wgs <dbl>,
## #   longitude_wgs <dbl>, notifikasi_status <fct>, race <fct>,
## #   kewarganegaraan <fct>, gender <fct>, nationality <fct>,
## #   klasifikasi_kes <fct>, cara_pengesanan_kes <fct>, jenis_rawatan <fct>,
## #   daerah <fct>, mukim <fct>, lokaliti <fct>, diagnosis2 <fct>, …

or, you can proceed with the analysis by removing the observation with no coordinates

kel_lepto2 <- kel_lepto %>% # save to different name to avoid losing the data
  filter(!is.na(latitude_wgs),
         !is.na(longitude_wgs))
glimpse(kel_lepto2)

## Rows: 1,106
## Columns: 26
## $ diagnosis                                  <fct> Leptospirosis, Leptospirosi…
## $ notifikasi_no                              <dbl> 2685725, 2728504, 2739963, …
## $ tahun_daftar                               <dbl> 2016, 2016, 2016, 2016, 201…
## $ epid_daftar                                <dbl> 6, 9, 11, 12, 18, 18, 19, 1…
## $ age                                        <dbl> 30, 23, 39, 43, 31, 34, 48,…
## $ alamat_semasa_kejadian                     <fct> "FELDA ARING", "LADANG U&I …
## $ poskod                                     <dbl> 18300, 18300, 18300, 18300,…
## $ latitude_rso                               <dbl> 478031, 459494, 441802, 488…
## $ longitude_rso                              <dbl> 548141, 564966, 547551, 547…
## $ latitude_wgs                               <dbl> 4.944824, 5.034372, 4.89743…
## $ longitude_wgs                              <dbl> 102.3491, 102.1477, 101.960…
## $ notifikasi_status                          <fct> Daftar Kes, Daftar Kes, Daf…
## $ race                                       <fct> Foreigner, Foreigner, Forei…
## $ kewarganegaraan                            <fct> Bukan Warganegara, Bukan Wa…
## $ gender                                     <fct> Male, Male, Male, Male, Mal…
## $ nationality                                <fct> INDONESIA, INDONESIA, INDON…
## $ klasifikasi_kes                            <fct> Sporadic, Sporadic, Sporadi…
## $ cara_pengesanan_kes                        <fct> Pasif, Pasif, Pasif, Pasif,…
## $ jenis_rawatan                              <fct> Wad Perubatan, Jabatan Kece…
## $ daerah                                     <fct> GUA MUSANG, GUA MUSANG, GUA…
## $ mukim                                      <fct> CHIKU, CHIKU, GALAS, CHIKU,…
## $ lokaliti                                   <fct> NA, "LADANG U & I, CIKU", "…
## $ diagnosis2                                 <fct> LEPTOSPIROSIS, LEPTOSPIROSI…
## $ sub_diagnosis                              <lgl> NA, NA, NA, NA, NA, NA, NA,…
## $ negeri                                     <fct> KELANTAN, KELANTAN, KELANTA…
## $ jenis_import_jangkitan_dalam_negara_negeri <fct> NA, NA, NA, NA, NA, NA, NA,…

Merge data

Merge population data to polygon

We want to match the two data based on MUKIM. Ensure the column containing the mukim is named the same for both data.

Use str() to look at the structure of the data.

str(kel_map)

## Classes 'sf' and 'data.frame':   66 obs. of  7 variables:
##  $ NEGERI    : chr  "KELANTAN" "KELANTAN" "KELANTAN" "KELANTAN" ...
##  $ DAERAH    : chr  "BACHOK" "BACHOK" "BACHOK" "BACHOK" ...
##  $ MUKIM     : chr  "BEKLAM" "GUNONG (GUNONG TIMOR)" "MAHLIGAI" "PERUPOK" ...
##  $ LELAKI    : int  4859 11100 4564 8777 9227 14140 5863 4929 18064 11645 ...
##  $ PEREMPUAN : num  4813 10884 4600 8614 8672 ...
##  $ JUM_JANTIN: num  9672 21984 9164 17391 17899 ...
##  $ geometry  :sfc_POLYGON of length 66; first list element: List of 1
##   ..$ : num [1:238, 1:2] 485502 485717 485966 486066 486068 ...
##   ..- attr(*, "class")= chr [1:3] "XY" "POLYGON" "sfg"
##  - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA
##   ..- attr(*, "names")= chr [1:6] "NEGERI" "DAERAH" "MUKIM" "LELAKI" ...

in ‘kel_map’ polygon data, the column name is ‘MUKIM’

str(kel_mukim)

## tibble [191 × 8] (S3: tbl_df/tbl/data.frame)
##  $ subdistrict: chr [1:191] "BEKLAM" "GUNONG (GUNONG TIMOR)" "MAHLIGAI" "PERUPOK" ...
##  $ 2016       : num [1:191] 11800 27100 11000 21100 23100 35400 17200 12500 45000 29800 ...
##  $ 2017       : num [1:191] 11800 27100 11000 21100 23100 35400 17200 12500 45000 29800 ...
##  $ 2018       : num [1:191] 12000 27500 11100 21400 23700 36200 17800 12800 45900 30500 ...
##  $ 2019       : num [1:191] 12200 28000 11300 21700 24200 37000 18400 13000 46800 31300 ...
##  $ 2020       : num [1:191] 12400 28400 11400 22000 24800 37800 19000 13300 47800 32000 ...
##  $ 2021       : num [1:191] 12400 28400 11400 22000 24800 37800 19000 13300 47800 32000 ...
##  $ 2022       : num [1:191] 12400 28400 11400 22000 24800 37800 19000 13300 47800 32000 ...

in ‘kel_mukim’ population data, the column name is ‘subdistrict’. We will rename the ‘subdistrict’ variable to ‘MUKIM’.

kel_mukim <- kel_mukim %>%
  rename(MUKIM = subdistrict)

str(kel_lepto2)

## tibble [1,106 × 26] (S3: tbl_df/tbl/data.frame)
##  $ diagnosis                                 : Factor w/ 1 level "Leptospirosis": 1 1 1 1 1 1 1 1 1 1 ...
##  $ notifikasi_no                             : num [1:1106] 2685725 2728504 2739963 2754040 2811959 ...
##  $ tahun_daftar                              : num [1:1106] 2016 2016 2016 2016 2016 ...
##  $ epid_daftar                               : num [1:1106] 6 9 11 12 18 18 19 17 21 23 ...
##  $ age                                       : num [1:1106] 30 23 39 43 31 34 48 55 48 39 ...
##  $ alamat_semasa_kejadian                    : Factor w/ 1073 levels "110, KESEDAR RENOK BARU",..: 108 689 1053 111 676 690 677 220 130 682 ...
##  $ poskod                                    : num [1:1106] 18300 18300 18300 18300 18300 ...
##  $ latitude_rso                              : num [1:1106] 478031 459494 441802 488502 496760 ...
##  $ longitude_rso                             : num [1:1106] 548141 564966 547551 547701 539072 ...
##  $ latitude_wgs                              : num [1:1106] 4.94 5.03 4.9 4.95 4.94 ...
##  $ longitude_wgs                             : num [1:1106] 102 102 102 102 102 ...
##  $ notifikasi_status                         : Factor w/ 2 levels "Abai Notifikasi",..: 2 2 2 2 2 2 2 2 2 2 ...
##  $ race                                      : Factor w/ 6 levels "Aborigines","Chinese",..: 3 3 3 3 3 3 3 3 3 3 ...
##  $ kewarganegaraan                           : Factor w/ 2 levels "Bukan Warganegara",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ gender                                    : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 1 2 1 ...
##  $ nationality                               : Factor w/ 8 levels "BANGLADESH","CAMBODIA",..: 4 4 4 4 3 3 3 8 4 6 ...
##  $ klasifikasi_kes                           : Factor w/ 2 levels "Outbreak","Sporadic": 2 2 2 2 2 2 2 2 2 2 ...
##  $ cara_pengesanan_kes                       : Factor w/ 2 levels "Aktif","Pasif": 2 2 2 2 2 2 2 2 2 2 ...
##  $ jenis_rawatan                             : Factor w/ 6 levels "ICU","Jabatan Kecemasan & Kemalangan",..: 6 2 6 6 6 6 6 6 6 6 ...
##  $ daerah                                    : Factor w/ 10 levels "BACHOK","GUA MUSANG",..: 2 2 2 2 2 2 2 10 2 2 ...
##  $ mukim                                     : Factor w/ 65 levels "ALOR PASIR","BADANG",..: 14 14 16 14 14 7 7 48 14 7 ...
##  $ lokaliti                                  : Factor w/ 710 levels "AIR BOL","AIR DERAS",..: NA 421 616 170 NA NA NA 667 323 NA ...
##  $ diagnosis2                                : Factor w/ 1 level "LEPTOSPIROSIS": 1 1 1 1 1 1 1 1 1 1 ...
##  $ sub_diagnosis                             : logi [1:1106] NA NA NA NA NA NA ...
##  $ negeri                                    : Factor w/ 1 level "KELANTAN": 1 1 1 1 1 1 1 1 1 1 ...
##  $ jenis_import_jangkitan_dalam_negara_negeri: Factor w/ 1 level "KELANTAN": NA NA NA NA NA NA NA NA NA NA ...

in ‘kel_lepto2’ population data, the column name is ‘mukim’. We will rename the ‘mukim’ variable to ‘MUKIM’.

kel_lepto2 <- kel_lepto2 %>%
  rename(MUKIM = mukim)

merge population data to polygon

kel <- merge(kel_map,kel_mukim, by.x="MUKIM", by.y="MUKIM", all.x=T, sort=F)
dim(kel)

## [1] 66 14

class(kel)

## [1] "sf"         "data.frame"

st_crs(kel)

## Coordinate Reference System:
##   User input: Kertau (RSO) / RSO Malaya (m) 
##   wkt:
## PROJCRS["Kertau (RSO) / RSO Malaya (m)",
##     BASEGEOGCRS["Kertau (RSO)",
##         DATUM["Kertau (RSO)",
##             ELLIPSOID["Everest 1830 (RSO 1969)",6377295.664,300.8017,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4751]],
##     CONVERSION["Rectified Skew Orthomorphic Malaya Grid (metres)",
##         METHOD["Hotine Oblique Mercator (variant A)",
##             ID["EPSG",9812]],
##         PARAMETER["Latitude of projection centre",4,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8811]],
##         PARAMETER["Longitude of projection centre",102.25,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8812]],
##         PARAMETER["Azimuth of initial line",323.0257905,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8813]],
##         PARAMETER["Angle from Rectified to Skew Grid",323.130102361111,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8814]],
##         PARAMETER["Scale factor on initial line",0.99984,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8815]],
##         PARAMETER["False easting",804670.24,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Engineering survey, topographic mapping."],
##         AREA["Malaysia - West Malaysia onshore."],
##         BBOX[1.21,99.59,6.72,104.6]],
##     ID["EPSG",3168]]

Check CRS

convert all point data to geometry (sf object). for WGS84, the CRS code is 4326.

lepto_wgs <- st_as_sf(kel_lepto2, 
                    coords = c("longitude_wgs", "latitude_wgs"), 
                    crs = 4326)

confirm CRS is WGS 84

st_crs(lepto_wgs)

## Coordinate Reference System:
##   User input: EPSG:4326 
##   wkt:
## GEOGCRS["WGS 84",
##     ENSEMBLE["World Geodetic System 1984 ensemble",
##         MEMBER["World Geodetic System 1984 (Transit)"],
##         MEMBER["World Geodetic System 1984 (G730)"],
##         MEMBER["World Geodetic System 1984 (G873)"],
##         MEMBER["World Geodetic System 1984 (G1150)"],
##         MEMBER["World Geodetic System 1984 (G1674)"],
##         MEMBER["World Geodetic System 1984 (G1762)"],
##         MEMBER["World Geodetic System 1984 (G2139)"],
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ENSEMBLEACCURACY[2.0]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     USAGE[
##         SCOPE["Horizontal component of 3D system."],
##         AREA["World."],
##         BBOX[-90,-180,90,180]],
##     ID["EPSG",4326]]

covert CRS to RSO (because the CRS for polygon map is in RSO-code 3168)

lepto_rso <- st_transform(lepto_wgs, 3168)

Check for outliers

plot the map to look for outliers

ggplot() +
  geom_sf(data = lepto_rso) +
  theme_bw()

if there are outliers, recheck the coordinates, update in the same linelisting and re-run the codes from beginning. or, you can proceed with analysis by excluding points outside the map’s boundaries.

select point only in Kelantan

lepto<- lepto_rso %>% 
  mutate(within_kel_map = lengths(st_within(lepto_rso, kel_map)))
lepto <- lepto %>% 
  filter(within_kel_map == 1)
glimpse(lepto)

## Rows: 1,106
## Columns: 26
## $ diagnosis                                  <fct> Leptospirosis, Leptospirosi…
## $ notifikasi_no                              <dbl> 2685725, 2728504, 2739963, …
## $ tahun_daftar                               <dbl> 2016, 2016, 2016, 2016, 201…
## $ epid_daftar                                <dbl> 6, 9, 11, 12, 18, 18, 19, 1…
## $ age                                        <dbl> 30, 23, 39, 43, 31, 34, 48,…
## $ alamat_semasa_kejadian                     <fct> "FELDA ARING", "LADANG U&I …
## $ poskod                                     <dbl> 18300, 18300, 18300, 18300,…
## $ latitude_rso                               <dbl> 478031, 459494, 441802, 488…
## $ longitude_rso                              <dbl> 548141, 564966, 547551, 547…
## $ notifikasi_status                          <fct> Daftar Kes, Daftar Kes, Daf…
## $ race                                       <fct> Foreigner, Foreigner, Forei…
## $ kewarganegaraan                            <fct> Bukan Warganegara, Bukan Wa…
## $ gender                                     <fct> Male, Male, Male, Male, Mal…
## $ nationality                                <fct> INDONESIA, INDONESIA, INDON…
## $ klasifikasi_kes                            <fct> Sporadic, Sporadic, Sporadi…
## $ cara_pengesanan_kes                        <fct> Pasif, Pasif, Pasif, Pasif,…
## $ jenis_rawatan                              <fct> Wad Perubatan, Jabatan Kece…
## $ daerah                                     <fct> GUA MUSANG, GUA MUSANG, GUA…
## $ MUKIM                                      <fct> CHIKU, CHIKU, GALAS, CHIKU,…
## $ lokaliti                                   <fct> NA, "LADANG U & I, CIKU", "…
## $ diagnosis2                                 <fct> LEPTOSPIROSIS, LEPTOSPIROSI…
## $ sub_diagnosis                              <lgl> NA, NA, NA, NA, NA, NA, NA,…
## $ negeri                                     <fct> KELANTAN, KELANTAN, KELANTA…
## $ jenis_import_jangkitan_dalam_negara_negeri <fct> NA, NA, NA, NA, NA, NA, NA,…
## $ geometry                                   <POINT [m]> POINT (484205.4 54689…
## $ within_kel_map                             <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, …

3. Disease Map

Map all the points on the Kelantan polygon map

ggplot() +
  geom_sf(data = kel_map) +
  geom_sf(data = lepto) +
    theme_bw()

alternatively, you can plot the points on open street view map

library(leaflet)
library(ggmap)

## ℹ Google's Terms of Service: <https://mapsplatform.google.com>
##   Stadia Maps' Terms of Service: <https://stadiamaps.com/terms-of-service/>
##   OpenStreetMap's Tile Usage Policy: <https://operations.osmfoundation.org/policies/tiles/>
## ℹ Please cite ggmap if you use it! Use `citation("ggmap")` for details.
## 
## Attaching package: 'ggmap'
## 
## 
## The following object is masked from 'package:plotly':
## 
##     wind

leaflet(data = lepto_wgs) %>%
  addTiles() %>%  # Add default OpenStreetMap map tiles
  addCircleMarkers(color = "red", radius = 1, fillOpacity = 0.7)

Map the disease by year.

Adjust the size of the plot using fig.width and fig.height. Adjust the size of the font by increasing or decreasing the text size.

facet_wrap() allows you to stratify the plot based on certain variables (eg: tahun_daftar). Change the variable name if you want to stratify based on other variables (eg: gender)

ggplot() +
  geom_sf(data = kel_map) +
  geom_sf(data = lepto) +
  ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022") +
  theme_bw() +
  facet_wrap(~ tahun_daftar) +
  theme(plot.title = element_text(size = 20),  strip.text = element_text(size = 20), axis.text.x=element_text(size=10), axis.text.y=element_text(size=10))

if you want to stratify based on more than one stratified variables (eg: gender and tahun_daftar), use facet_grid()

ggplot() +
  geom_sf(data = kel_map) +
  geom_sf(data = lepto) +
  ggtitle("Map of Leptospirosis Cases in Kelantan by gender for 2016-2022") +
  theme_bw() +
  facet_grid(gender ~ tahun_daftar) +
  theme(plot.title = element_text(size = 24),  strip.text = element_text(size = 20), axis.text.x=element_text(size=10), axis.text.y=element_text(size=10))

Interactive plot using plotly

library(plotly)

# Example ggplot object without facet_wrap
p <- ggplot() +
  geom_sf(data = kel_map) +
  geom_sf(data = lepto) +
  ggtitle("Map of Leptospirosis Cases in Kelantan for 2016-2022") +
  theme_bw() +
  theme(
    plot.title = element_text(size = 20),
    axis.text.x = element_text(size = 10),
    axis.text.y = element_text(size = 10)
  )

# Convert the ggplot object to plotly
ggplotly(p)

Interactive map according to year using shiny

library(shiny)

# Define UI for the application
ui <- fluidPage(
  titlePanel("Interactive Map of Leptospirosis Cases in Kelantan"),
  
  # Sidebar with two dropdowns: one for year, one for daerah
  sidebarLayout(
    sidebarPanel(
      # Dropdown to select year
      selectInput("year", "Select Year:", choices = unique(lepto$tahun_daftar)),
      
      # Dropdown to select daerah
      selectInput("daerah", "Select Daerah:", choices = unique(lepto$daerah))
    ),
    
    # Display the interactive plot
    mainPanel(
      plotlyOutput("leptoPlot")
    )
  )
)

# Define server logic
server <- function(input, output) {
  output$leptoPlot <- renderPlotly({
    
    # Filter data based on the selected year and daerah
    filtered_lepto <- lepto[lepto$tahun_daftar == input$year & lepto$daerah == input$daerah, ]
    
    # Create the ggplot
    p <- ggplot() +
      geom_sf(data = kel_map) +
      geom_sf(data = filtered_lepto) +
      ggtitle(paste("Map of Leptospirosis Cases in", input$daerah, "for", input$year)) +
      theme_bw() +
      theme(
        plot.title = element_text(size = 20),
        axis.text.x = element_text(size = 10),
        axis.text.y = element_text(size = 10)
      )
    
    # Convert ggplot to plotly
    ggplotly(p)
  })
}

# Run the application
shinyApp(ui = ui, server = server)

4. Point Pattern Analysis

The points need to be converted to spatial data

As we are goint to analyze the disease for each year, lets extract the leptospirosis cases by year

lepto16 <- subset(lepto, tahun_daftar=="2016")
lepto17 <- subset(lepto, tahun_daftar=="2017")
lepto18 <- subset(lepto, tahun_daftar=="2018")
lepto19 <- subset(lepto, tahun_daftar=="2019")
lepto20 <- subset(lepto, tahun_daftar=="2020")
lepto21 <- subset(lepto, tahun_daftar=="2021")
lepto22 <- subset(lepto, tahun_daftar=="2022")

set observation window

kel_map.owin <- as.owin(kel_map)
plot(kel_map.owin)

Convert to spatial object

lepto16.sp <- as(lepto16, "Spatial")
lepto17.sp <- as(lepto17, "Spatial")
lepto18.sp <- as(lepto18, "Spatial")
lepto19.sp <- as(lepto19, "Spatial")
lepto20.sp <- as(lepto20, "Spatial")
lepto21.sp <- as(lepto21, "Spatial")
lepto22.sp <- as(lepto22, "Spatial")

Convert to planar point pattern (ppp) object

coords16 <- coordinates(lepto16.sp)  
lepto16.ppp <- ppp(x = coords16[,1], y = coords16[,2], window = kel_map.owin)

## Warning: data contain duplicated points

coords17 <- coordinates(lepto17.sp)  
lepto17.ppp <- ppp(x = coords17[,1], y = coords17[,2], window = kel_map.owin)

coords18 <- coordinates(lepto18.sp)  
lepto18.ppp <- ppp(x = coords18[,1], y = coords18[,2], window = kel_map.owin)

## Warning: data contain duplicated points

coords19 <- coordinates(lepto19.sp)  
lepto19.ppp <- ppp(x = coords19[,1], y = coords19[,2], window = kel_map.owin)

coords20 <- coordinates(lepto20.sp)  
lepto20.ppp <- ppp(x = coords20[,1], y = coords20[,2], window = kel_map.owin)

coords21 <- coordinates(lepto21.sp)  
lepto21.ppp <- ppp(x = coords21[,1], y = coords21[,2], window = kel_map.owin)

coords22 <- coordinates(lepto22.sp)  
lepto22.ppp <- ppp(x = coords22[,1], y = coords22[,2], window = kel_map.owin)

## Warning: data contain duplicated points

4.1 Density-based Analysis

Quadrat analysis

Quadrat count

par( mfrow= c(2,4) ) #combine all plot in 1 view, 2 row, 4 column
# 2016
quadr_count_lepto16 <- quadratcount(lepto16.ppp, 
                                  nx = 10,
                                  ny = 14)
plot(lepto16.ppp, pch = 20, cex = 0.1, main = "2016")
plot(quadr_count_lepto16, add = TRUE, cex = 2)

# 2017
quadr_count_lepto17 <- quadratcount(lepto17.ppp, 
                                  nx = 10,
                                  ny = 14)
plot(lepto17.ppp, pch = 20, cex = 0.1, main = "2017")
plot(quadr_count_lepto17, add = TRUE, cex = 2)

# 2018
quadr_count_lepto18 <- quadratcount(lepto18.ppp, 
                                  nx = 10,
                                  ny = 14)
plot(lepto18.ppp, pch = 20, cex = 0.1, main = "2018")
plot(quadr_count_lepto18, add = TRUE, cex = 2)

# 2019
quadr_count_lepto19 <- quadratcount(lepto19.ppp, 
                                  nx = 10,
                                  ny = 14)
plot(lepto19.ppp, pch = 20, cex = 0.1, main = "2019")
plot(quadr_count_lepto19, add = TRUE, cex = 2)

# 2020
quadr_count_lepto20 <- quadratcount(lepto20.ppp, 
                                  nx = 10,
                                  ny = 14)
plot(lepto20.ppp, pch = 20, cex = 0.1, main = "2020")
plot(quadr_count_lepto20, add = TRUE, cex = 2)

# 2021
quadr_count_lepto21 <- quadratcount(lepto21.ppp, 
                                  nx = 10,
                                  ny = 14)
plot(lepto21.ppp, pch = 20, cex = 0.1, main = "2021")
plot(quadr_count_lepto21, add = TRUE, cex = 2)

# 2022
quadr_count_lepto22 <- quadratcount(lepto22.ppp, 
                                  nx = 10,
                                  ny = 14)
plot(lepto22.ppp, pch = 20, cex = 0.1, main = "2022")
plot(quadr_count_lepto22, add = TRUE, cex = 2)

Test for CSR

Chi-squared test

Chi-squared goodness-of-fit test

chi_lepto16 <- quadrat.test(lepto16.ppp, nx= 10, ny=14)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

chi_lepto17 <- quadrat.test(lepto17.ppp, nx= 10, ny=14)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

chi_lepto18 <- quadrat.test(lepto18.ppp, nx= 10, ny=14)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

chi_lepto19 <- quadrat.test(lepto19.ppp, nx= 10, ny=14)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

chi_lepto20 <- quadrat.test(lepto20.ppp, nx= 10, ny=14)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

chi_lepto21 <- quadrat.test(lepto21.ppp, nx= 10, ny=14)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

chi_lepto22 <- quadrat.test(lepto22.ppp, nx= 10, ny=14)

## Warning: Some expected counts are small; chi^2 approximation may be inaccurate

Display result

chi_lepto16.df <- data.frame(Dataset = "lepto16.ppp",
                       TestStatistic = chi_lepto16$statistic,
                       PValue = chi_lepto16$p.value)
chi_lepto17.df <- data.frame(Dataset = "lepto17.ppp",
                       TestStatistic = chi_lepto17$statistic,
                       PValue = chi_lepto17$p.value)
chi_lepto18.df <- data.frame(Dataset = "lepto18.ppp",
                       TestStatistic = chi_lepto18$statistic,
                       PValue = chi_lepto18$p.value)
chi_lepto19.df <- data.frame(Dataset = "lepto19.ppp",
                       TestStatistic = chi_lepto19$statistic,
                       PValue = chi_lepto19$p.value)
chi_lepto20.df <- data.frame(Dataset = "lepto20.ppp",
                       TestStatistic = chi_lepto20$statistic,
                       PValue = chi_lepto20$p.value)
chi_lepto21.df <- data.frame(Dataset = "lepto21.ppp",
                       TestStatistic = chi_lepto21$statistic,
                       PValue = chi_lepto21$p.value)
chi_lepto22.df <- data.frame(Dataset = "lepto22.ppp",
                       TestStatistic = chi_lepto22$statistic,
                       PValue = chi_lepto22$p.value)
chi_quadlepto <- bind_rows(
  data.frame(Dataset = "2016", chi_lepto16.df), 
  data.frame(Dataset = "2017", chi_lepto17.df), 
  data.frame(Dataset = "2018", chi_lepto18.df), 
  data.frame(Dataset = "2019", chi_lepto19.df), 
  data.frame(Dataset = "2020", chi_lepto20.df), 
  data.frame(Dataset = "2021", chi_lepto21.df), 
  data.frame(Dataset = "2022", chi_lepto22.df))
  chi_quadlepto

##        Dataset   Dataset.1 TestStatistic        PValue
## X2...1    2016 lepto16.ppp     1171.8180 6.181358e-181
## X2...2    2017 lepto17.ppp      463.4975  3.361411e-47
## X2...3    2018 lepto18.ppp      479.9587  5.079450e-50
## X2...4    2019 lepto19.ppp      312.4394  6.684179e-23
## X2...5    2020 lepto20.ppp      202.7419  2.378868e-08
## X2...6    2021 lepto21.ppp      173.7245  2.476629e-05
## X2...7    2022 lepto22.ppp      446.4550  2.619867e-44

Monte Carlo based test

Monte Carlo test if chi-squared test assumption not met

# run monte carlo test

mc_qlepto16 <- quadrat.test(lepto16.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto17 <- quadrat.test(lepto17.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto18 <- quadrat.test(lepto18.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto19 <- quadrat.test(lepto19.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto20 <- quadrat.test(lepto20.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto21 <- quadrat.test(lepto21.ppp, nx= 10, ny=14, method = "MonteCarlo")
mc_qlepto22 <- quadrat.test(lepto22.ppp, nx= 10, ny=14, method = "MonteCarlo")

# convert to data frame
mc_qlepto16.df <- data.frame(Dataset = "lepto16.ppp",
                       TestStatistic = mc_qlepto16$statistic,
                       PValue = mc_qlepto16$p.value)
mc_qlepto17.df <- data.frame(Dataset = "lepto17.ppp",
                       TestStatistic = mc_qlepto17$statistic,
                       PValue = mc_qlepto17$p.value)
mc_qlepto18.df <- data.frame(Dataset = "lepto18.ppp",
                       TestStatistic = mc_qlepto18$statistic,
                       PValue = mc_qlepto18$p.value)
mc_qlepto19.df <- data.frame(Dataset = "lepto19.ppp",
                       TestStatistic = mc_qlepto19$statistic,
                       PValue = mc_qlepto19$p.value)
mc_qlepto20.df <- data.frame(Dataset = "lepto20.ppp",
                       TestStatistic = mc_qlepto20$statistic,
                       PValue = mc_qlepto20$p.value)
mc_qlepto21.df <- data.frame(Dataset = "lepto21.ppp",
                       TestStatistic = mc_qlepto21$statistic,
                       PValue = mc_qlepto21$p.value)
mc_qlepto22.df <- data.frame(Dataset = "lepto22.ppp",
                       TestStatistic = mc_qlepto22$statistic,
                       PValue = mc_qlepto22$p.value)

# combine rows
mc_qlepto <- bind_rows(
  data.frame(Dataset = "2016", mc_qlepto16.df), 
  data.frame(Dataset = "2017", mc_qlepto17.df), 
  data.frame(Dataset = "2018", mc_qlepto18.df), 
  data.frame(Dataset = "2019", mc_qlepto19.df), 
  data.frame(Dataset = "2020", mc_qlepto20.df), 
  data.frame(Dataset = "2021", mc_qlepto21.df), 
  data.frame(Dataset = "2022", mc_qlepto22.df))
mc_qlepto

##        Dataset   Dataset.1 TestStatistic PValue
## X2...1    2016 lepto16.ppp     1171.8180  0.002
## X2...2    2017 lepto17.ppp      463.4975  0.003
## X2...3    2018 lepto18.ppp      479.9587  0.004
## X2...4    2019 lepto19.ppp      312.4394  0.005
## X2...5    2020 lepto20.ppp      202.7419  0.010
## X2...6    2021 lepto21.ppp      173.7245  0.033
## X2...7    2022 lepto22.ppp      446.4550  0.007

A significant (p-value<0.05) chi-square or Monte-Carlo test result would indicate that the points are clustered and non-randomly distributed, suggesting the presence of spatial processes such as spatial contagion, spatial dependence, or spatial interaction.

Intensity analysis

par( mfrow= c(2,4) )

inten_lepto16 <-intensity(quadr_count_lepto16)
plot(intensity(quadr_count_lepto16, image = TRUE), main = "2016", las = 1)
plot(lepto16, pch = 20, cex = 0.6, add = TRUE)

## Warning in plot.sf(lepto16, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute

inten_lepto17 <-intensity(quadr_count_lepto17)
plot(intensity(quadr_count_lepto17, image = TRUE), main = "2017", las = 1)
plot(lepto17, pch = 20, cex = 0.6, add = TRUE)

## Warning in plot.sf(lepto17, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute

inten_lepto18 <-intensity(quadr_count_lepto18)
plot(intensity(quadr_count_lepto18, image = TRUE), main = "2018", las = 1)
plot(lepto18, pch = 20, cex = 0.6, add = TRUE)

## Warning in plot.sf(lepto18, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute

inten_lepto19 <-intensity(quadr_count_lepto19)
plot(intensity(quadr_count_lepto19, image = TRUE), main = "2019", las = 1)
plot(lepto19, pch = 20, cex = 0.6, add = TRUE)

## Warning in plot.sf(lepto19, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute

inten_lepto20 <-intensity(quadr_count_lepto20)
plot(intensity(quadr_count_lepto20, image = TRUE), main = "2020", las = 1)
plot(lepto20, pch = 20, cex = 0.6, add = TRUE)

## Warning in plot.sf(lepto20, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute

inten_lepto21 <-intensity(quadr_count_lepto21)
plot(intensity(quadr_count_lepto21, image = TRUE), main = "2021", las = 1)
plot(lepto21, pch = 20, cex = 0.6, add = TRUE)

## Warning in plot.sf(lepto21, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute

inten_lepto22 <-intensity(quadr_count_lepto22)
plot(intensity(quadr_count_lepto22, image = TRUE), main = "2022", las = 1)
plot(lepto22, pch = 20, cex = 0.6, add = TRUE)

## Warning in plot.sf(lepto22, pch = 20, cex = 0.6, add = TRUE): ignoring all but
## the first attribute

mtext("Intensity Maps of Leptospirosis Cases in Kelantan 2016-2022", side = 1, line = -1, cex = 1.5, outer = TRUE)

###Kernel Density Estimation (KDE)

pre-determined bandwidth (e.g. 5000metres)

kde.lepto16 <- density(lepto16.ppp, sigma = 5000) #2016
kde.lepto17 <- density(lepto17.ppp, sigma = 5000) #2017
kde.lepto18 <- density(lepto18.ppp, sigma = 5000) #2018
kde.lepto19 <- density(lepto19.ppp, sigma = 5000) #2019
kde.lepto20 <- density(lepto20.ppp, sigma = 5000) #2020
kde.lepto21 <- density(lepto21.ppp, sigma = 5000) #2021
kde.lepto22 <- density(lepto22.ppp, sigma = 5000) #2022

par( mfrow= c(2,4) )
plot(kde.lepto16, main = "2016", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto16, add = TRUE)
plot(kde.lepto17, main = "2017", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto17, add = TRUE)
plot(kde.lepto18, main = "2018", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto18, add = TRUE)
plot(kde.lepto19, main = "2019", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto19, add = TRUE)
plot(kde.lepto20, main = "2020", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto20, add = TRUE)
plot(kde.lepto21, main = "2021", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto21, add = TRUE)
plot(kde.lepto22, main = "2022", cex.main = 1.5, font.main = 2, las = 1)
contour(kde.lepto22, add = TRUE)

mtext("Kernel Density Estimate (KDE) Heatmaps of Leptospirosis Cases in Kelantan 2016-2022", side = 1, line = -1, cex = 2, outer = TRUE)

automated bandwidth KDE using bw selector Likelihood Cross Validation

kde.lepto16bw <- density(lepto16.ppp, sigma = bw.ppl(lepto16.ppp))
kde.lepto17bw <- density(lepto17.ppp, sigma = bw.ppl(lepto17.ppp))
kde.lepto18bw <- density(lepto18.ppp, sigma = bw.ppl(lepto18.ppp))
kde.lepto19bw <- density(lepto19.ppp, sigma = bw.ppl(lepto19.ppp))
kde.lepto20bw <- density(lepto20.ppp, sigma = bw.ppl(lepto20.ppp))
kde.lepto21bw <- density(lepto21.ppp, sigma = bw.ppl(lepto21.ppp))
kde.lepto22bw <- density(lepto22.ppp, sigma = bw.ppl(lepto22.ppp))

par( mfrow= c(2,4) )
plot(kde.lepto16bw, main = "2016", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto16bw, add = TRUE)
plot(kde.lepto17bw, main = "2017", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto17bw, add = TRUE)
plot(kde.lepto18bw, main = "2018", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto18bw, add = TRUE)
plot(kde.lepto19bw, main = "2019", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto19bw, add = TRUE)
plot(kde.lepto20bw, main = "2020", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto20bw, add = TRUE)
plot(kde.lepto21bw, main = "2021", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto21bw, add = TRUE)
plot(kde.lepto22bw, main = "2022", cex.main = 1.5, cex.axis = 2, font.main = 2, las = 1)
contour(kde.lepto22bw, add = TRUE)

mtext("Kernel Density Estimate (KDE) Heatmaps of Leptospirosis Cases in Kelantan 2016-2022", side = 1, line = -1, cex = 2, outer = TRUE)

Peaks: Peaks in the KDE plot represent areas where data points are concentrated, indicating higher density or probability of observations in those regions. The height of a peak shows how likely it is to find data points within that range.

Spread/Width: The width or spread of the KDE curve indicates the variability of the data. A wider curve suggests that the data points are more spread out, while a narrower curve indicates that the data points are concentrated in a smaller range.

Shape: The shape of the KDE plot can give insights into the distribution of your data. For instance:

Unimodal: A single peak suggests the data follows a normal-like distribution or another unimodal distribution.

Bimodal/Multimodal: Multiple peaks indicate that the data may have more than one mode, or there are distinct groups or clusters in the dataset.

The lighter the color, the higher the density

4.2 Distance-based Analysis

Average Nearest Neighbor (ANN) Analysis

The Average Nearest Neighbor (ANN) method measures the distance from each point in a point pattern to its nearest neighboring point and calculates the average of these distances.

This method focuses on local point-to-point relationships and is used to detect clustering or dispersion at a single scale (i.e., based on the nearest neighbor distance).

ann16 <- round(mean(nndist(lepto16.ppp)), 2)
ann17 <- round(mean(nndist(lepto17.ppp)), 2)
ann18 <- round(mean(nndist(lepto18.ppp)), 2)
ann19 <- round(mean(nndist(lepto19.ppp)), 2)
ann20 <- round(mean(nndist(lepto20.ppp)), 2)
ann21 <- round(mean(nndist(lepto21.ppp)), 2)
ann22 <- round(mean(nndist(lepto22.ppp)), 2)

annual_ann <- data.frame(
  Year = c(2016, 2017, 2018, 2019, 2020, 2021, 2022),
  Nearest_Neighbor_Distance = c(ann16, ann17, ann18, ann19, ann20, ann21, ann22)
)
print(annual_ann)

##   Year Nearest_Neighbor_Distance
## 1 2016                   1793.08
## 2 2017                   3282.28
## 3 2018                   2792.46
## 4 2019                   4385.30
## 5 2020                   7358.82
## 6 2021                   6394.20
## 7 2022                   2533.81

Nearest Neigbour Index (NNI)

If the index (average nearest neighbor ratio) is less than 1, the pattern exhibits clustering. If the index is greater than 1, the trend is toward dispersion.

nni16 <-nni(lepto16)

## Warning: data contain duplicated points

nni17 <-nni(lepto17)
nni18 <-nni(lepto18)

## Warning: data contain duplicated points

nni19 <-nni(lepto19)
nni20 <-nni(lepto20)
nni21 <-nni(lepto21)
nni22 <-nni(lepto22)

## Warning: data contain duplicated points

annual_nni <- bind_rows(
  data.frame(Dataset = "2016", nni16),
  data.frame(Dataset = "2017", nni17),
  data.frame(Dataset = "2018", nni18),
  data.frame(Dataset = "2019", nni19),
  data.frame(Dataset = "2020", nni20),
  data.frame(Dataset = "2021", nni21),
  data.frame(Dataset = "2022", nni22)
)
annual_nni

##   Dataset       NNI     z.score            p expected.mean.distance
## 1    2016 0.5974538 -14.6318827 1.758528e-48               3001.196
## 2    2017 0.8069370  -4.6425708 3.441007e-06               4067.577
## 3    2018 0.6811086  -8.1163418 4.804463e-16               4099.872
## 4    2019 0.7605712  -4.5112102 6.445880e-06               5765.793
## 5    2020 1.0425996   0.5281551 5.973917e-01               7058.142
## 6    2021 0.9495579  -0.5948619 5.519358e-01               6733.872
## 7    2022 0.6918232  -8.9993010 2.271594e-19               3662.505
##   observed.mean.distance
## 1               1793.076
## 2               3282.278
## 3               2792.458
## 4               4385.297
## 5               7358.816
## 6               6394.201
## 7               2533.806

Functions

G-function

Empirical values greater than theoretical (Poisson) values suggest clustering. Envelope denote simulations to test for CSR

par( mfrow= c(2,4) )
G_lepto16 <- plot(envelope(lepto16.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto17 <- plot(envelope(lepto17.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto18 <- plot(envelope(lepto18.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto19 <- plot(envelope(lepto19.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto20 <- plot(envelope(lepto20.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto21 <- plot(envelope(lepto21.ppp, Gest, nsim = 99, verbose = FALSE))
G_lepto22 <- plot(envelope(lepto22.ppp, Gest, nsim = 99, verbose = FALSE))

K-function

Empirical values greater than theoretical (Poisson) values suggest clustering. Envelope denote simulations to test for CSR Takes long time, for demonstration purpose, we will set the Monte Carlo simulation to only 3

par( mfrow= c(2,4) )
K_lepto16 <- plot(envelope(lepto16.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto17 <- plot(envelope(lepto17.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto18 <- plot(envelope(lepto18.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto19 <- plot(envelope(lepto19.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto20 <- plot(envelope(lepto20.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto21 <- plot(envelope(lepto21.ppp, Kest, nsim = 3, verbose = FALSE))
K_lepto22 <- plot(envelope(lepto22.ppp, Kest, nsim = 3, verbose = FALSE))

other summary function - Lest, Fest, pcf

Spatial Workshop

Hazlienor

2024-09-10